Field-induced vacancy localization in a driven lattice gas: Scaling of steady states. 
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With the help of Monte Carlo simulations and a mean-field theory, we investigate the ordered 
{3J[)' steady-state structures resulting from the motion of a single vacancy on a periodic lattice which is 

^ , filled with two species of oppositely "charged" particles. An external field biases particle-vacancy 

■ exchanges according to the particle's charge, subject to an excluded volume constraint. The steady 

state exhibits charge segregation, and the vacancy is localized at one of the two characteristic 
interfaces. Charge and hole density profiles, an appropriate order parameter and the interfacial 
regions themselves exhibit characteristic scaling properties with system size and field strength. The 
lattice spacing is found to play a significant role within the mean-field theory. 



I. INTRODUCTION 



Systems in non-equilibrium steady states have attracted considerable interest in the past decade While pre- 
senting a wealth of unexpected, intriguing phenomena, they are still quite poorly understood at a fundamental level. 
^ , It is therefore natural to investigate simple model systems to identify generic behavior, before turning to real systems 
' which are usually far more complex. 

a A particularly interesting class of model systems is based on lattice gas models, involving one or several species of 
particles whose motion is biased in a specified direction. If the boundary conditions are open or periodic, the bias 
can drive the system out of a well-known equilibrium state, into novel non-equilibrium steady states which typically 
(~| . carry global particle currents. Characteristic configurations, particle-particle correlations and even phase transitions 
O tend to be profoundly affected by the bias. Equilibrium phases can be suppressed, universality classes may change 
, O, , and entirely new transitions can emerge. For example, in a simple driven Ising lattice gas with periodic boundary 
conditions , the bias suppresses one of the two ground states of the equilibrium system and fundamentally changes 
T-H , the universal properties of the Ising order-disorder transition [^-||. In the high-temperature phase, it induces generic 
K*' ' long-range correlations 0, which characterize all models of this type. Other anomalies are observed below Tc 

' If the Ising symmetry is generalized [^,0 to include two (or more) species of particles which respond differently 
to the drive, such systems will generically exhibit blocking transitions, similar to traffic jams, in which one species 
impedes the motion of the other. These instabilities are genuine non-equilibrium transitions: they do not exist in 
. the equilibrium limit and are controlled by particle density and bias strength, rather than temperature. The ordered 
' phases exhibit characteristic spatial structures. Related real and model systems include water-in-oil microcmulsions 
^\ , in external electric fields , gel electrophoresis [|l3| and traffic flow . 

' In this paper, we focus on a three-state lattice gas consisting of holes and two distinct species of particles driven 
I by an external field E (the bias) in opposite directions . We name the two species "positive" and "negative" , in 
H • analogy to charged particles in an electric field. The bias clearly breaks the Potts symmetry of the stochastic variable, 
by acting differently on each species. The only interaction between the particles is an excluded volume constraint 
so that (i) each site can be occupied by at most one particle, and (ii) particle-particle ("charge") exchanges are not 
allowed. In the absence of other interparticle interactions, the temperature dependence of the system, reflecting a 
Q . coupling to a heat bath, may be absorbed into the drive. Hence, the model is a high-field, high-temperature limit of 
a more complicated interacting system. 

On a fully periodic square lattice, this system undergoes a blocking transition controlled by field strength and 
particle density, separating a homogeneous phase from a spatially inhomogeneous one [ p5| . For small mass density 
5h ' and drive, the steady-state configurations are disordered, so that both particle densities are homogeneous and a 
significant charge current persists. In contrast, if a threshold mass density is exceeded, the particles form a single 
compact strip transverse to the field while the rest of the lattice remains essentially empty. The particle-rich region 
itself consists of two strips, also oriented transverse to the field, each dominated by one single species. In this phase, 
the particles impede one another, due to the excluded volume constraint, so that the charge current is much smaller. 
Other ordered phases, with nonzero winding number around the lattice ("barber poles"), are observed in systems 
with rectangular aspect ratios . An analytical solution in the frame of a mean-field theory was presented in 
Ref. . With a slight geometrical modification, the model was also investigated by Foster and Godreche . 
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Here, we focus on a novel aspect of the blocking transition, namely, a localization phenomenon occurring in systems 
near complete filling ||l9t| . Thus, all lattice sites except a single vacant one are occupied by particles. For simplicity, 
we consider the symmetric situation, i.e., the particle numbers of each species differ at most by one. Starting from a 
disordered initial configuration, particles may exchange only with the vacancy. As a result, the hole diffuses through 
the lattice. However, it does not perform a Brownian random walk, since the jump rate for a particle-hole exchange 
depends on the charge and direction of motion of the particle. By virtue of the bias, positive and negative particles 
are transported in opposite directions: The two particle species eventually segregate, provided the field exceeds a 
certain threshold, corresponding to the transition line |p^ . When the steady state is reached, two strips have formed, 
filled by positive and negative particles, respectively. The hole itself ends up "trapped" on one of the two interfaces 
between the two ordered regions. Its location is the remnant of the empty region observed at finite hole density. 

This problem, in both its static and dynamic aspects, is an example of a much-wider ranging class of interacting 
random walk and defect-mediated domain growth problems. The hole is a random walker whose motion changes its 
environment, but the environment reacts by determining the local jump rates. The vacancy plays the role of a highly 
mobile defect poj , interacting with an otherwise immobile background. The time evolution of the system, from an 
initially disordered particle background to two ordered strips, poses a domain growth problem [pl| . Clearly, a good 
understanding of the final steady states and their associated scaling properties is the first step in the analysis of the 
ordering process. This study forms the subject of this paper. We report elsewhere on the full dynamics ||22| . 

The key results of our study p^ ] can be summarized as follows. First, we establish and confirm the characteristic 
scaling forms of the order parameter and the density profiles. Further, focusing only on the interfacial (as opposed to 
the fully ordered) regions of the profiles, we find that both interfaces are independent of the longitudinal system size 
and that their widths are controlled by the drive alone. These findings are reflected in our mean-field theory. Our 
results are limited in two ways: first, by the onset of the phase transition for small E, and second, by the breakdown 
of the naive continuum limit for large E. 

This paper is organized as follows. In Section 2, we give a precise definition of the microscopic model which underlies 
the Monte Carlo simulations. The relevant control and order parameters are defined. To set the scene, we provide 
a brief summary of earlier work. In particular, we discuss the blocking transition and its description in terms of a 
mean- field theory. In Section 3, we investigate the scaling properties of the order parameter and the profiles, based 
on Monte Carlo simulations and the exact solution of the mean-field equations. We conclude with a summary and 
some comments. 



II. THE MODEL: MICROSCOPICS AND MEAN-FIELD THEORY 



In this section, we provide the necessary background. We begin with the microscopic definition of the model, 
followed by a summary of its phenomenology. We then provide a different perspective, by sketching the mean-field 
theory and its main results. We close with some technical details of the simulations. 

Our model is defined on a two-dimensional square lattice of x Ly sites with fully periodic boundary conditions. 
Each site, except one, can be occupied by a positive or negative particle. The remaining site is left empty. The 
resulting configurations can be described by an occupation variable n+y (n^y), taking the value +1 if a positive 
(negative) charge is present at site (x, y) and zero otherwise. This enforces the excluded volume constraint. There are 
no other interactions between the particles. Turning to dynamics, particles may jump only onto the vacant site. In 
the absence of the external field, the vacancy exchanges randomly with any of its four nearest neighbors, independent 
of their charge and the direction of the move. This symmetry is broken by the "electric" field E, which is chosen to 
be uniform in space and time and directed along the positive y-axis. For nonzero E, jumps transverse to the field are 
still random; however, parallel jumps are now biased: positive (negative) charges jump preferentially along (against) 
E. Specifically, the exchange rate of the hole with a randomly chosen nearest neighbor is given by the Metropolis 
rate Q: 

W = T[im{l,exp{qESy)} (1) 

where q — +1 (—1) for a positive (negative) particle and Sy = 0, ±a is the change of the y-coordinate of the particle 
due to the jump. This choice mimics the local energetics of charges in a uniform field. The lattice constant a will be 
set to 1. 

The dynamics of the model can be summarized by a master equation p5| , for the probability P(C,t) to find the 
system in the configuration C = n~y^ at time t: 

^P(C, t)=J2 {WiC ^ C)P(C', t) - W{C ^ C')PiC, t)} . (2) 
c 
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Here, W{C — > C) is the transition rate from C to C", specified by Eqn For E < oo, P{C,t) approaches a unique 
steady-state solution P*{C) in the hmit t ^ oo. For closed boundary conditions, P*{C) follows from equilibrium 
statistical mechanics, being the Boltzmann factor of a system of noninteracting charges in a uniform field. For periodic 
boundary conditions, however, there is no uniquely defined static potential for i? 7^ so that P*{C) is not a priori 
known. Instead, it has to be found from an explicit solution of Eqn (Q). Unfortunately, such solutions are available 
only for a few, mostly one-dimensional, cases. Here, we only know the E ^ solution: the system is again in 
equilibrium, the particles diffuse randomly and P*{C) is independent of configuration, i.e., P* cx 1. 

The control parameters of this model are easily identified. In addition to the driving field E and the system size, 
Lx y. Ly, we can adjust the mass density 
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LxLy 



nty + nxy) (3) 



as well as the net charge density of the system 

P 



LxLy 



Since the particle number of each species is separately conserved, both densities are also conserved. For our case, 
there is always a single hole, so that the particle density 

m = 1 - -i- (5) 

J^xJ^y 

depends on the system size. Since we focus on nearly equal numbers of positive and negative particles, the net charge 
density is zero for systems with an odd number of sites and —\/ [LxLy) for an even number (the hole always takes 
the place of a positive particle). In the simulations, this small difference does not appear to lead to observable effects, 
unlike the case oi p = 0{1) Pq| . 

A brief description of the blocking transition and the associated phases will be helpful. For small values of drive 
and total mass, the system is in the disordered phase, characterized by spatially uniform mass and charge densities. 
A significant charge current flows in this phase. As or m increase, a transition into an ordered phase, with 
spatially inhomogeneous densities, occurs. For systems with aspect ratios near unity, each species of particles forms 
a compact, stable strip transverse to the electric field. The strip of positive charges is located directly "upfield" from 
the negative strip, so that the strips block each other, due to the excluded volume constraint. The rest of the lattice 
remains essentially empty. Clearly, the current is much smaller in this phase. In the following, we will investigate the 
structure of these transverse strips when the empty region has shrunk to a single hole. We never observe strips with 
nonzero winding number: they appear to be suppressed near complete filling. 

To distinguish ordered and disordered phases, a suitable order parameter is needed. It is convenient to introduce 
the local hole and charge densities: 

(t>x,y = l- {nt^y + n-y) and V'x.j, = n+j, - n'^. (6) 

Since our system does not develop inhomogeneities in the x-direction, it is natural to focus on the mass and charge 
density profiles: 

Hy) ^ ^Y<l>x,y and ip{y) = ^^'ipx,y (7) 

^ X ^ X 

Following Rcf. we define an order parameter 
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The angular bracket denotes a configurational average. Squaring the charge density profile (which can have either 
sign) prevents unwanted cancellations in the sum over y. In the ordered phase, Ql is 0(1), while being only order 
0{l/{mLx)) in the disordered phase. Roughly speaking, mLyQi^ counts the ordered rows transverse to the external 
field. For a perfectly ordered system, Ql would be unity. Clearly, other definitions of an order parameter are possible. 
In particular, the amplitude of the lowest Fourier component of either tp{y) or (j){y) is a much more sensitive measure 
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for a study of the transition line |^JT^. Here, however, our focus is not on the transition, but on the structure of 
ordered states, so that Ql serves its purpose well. 

Finally, let us add a comment on the transition line. Earlier simulation data show that the threshold mass, 
mc, depends str ongly o n E and the longitudinal system length Ly, but only weakly (if at all) on L^. It can be first 
or second order |17| , p9[ , in different regions of parameter space. In our case, where only a single hole is present, the 
mass density is unity, to excellent accuracy. Therefore, the transition is controlled by E and Ly alone. 

We now turn to a brief summary of the theoretical analysis ||l^,|l^ which will be essential for the following. Even 
though the master equation is just a linear equation, in practice it is not susceptible to theoretical analysis. To 
proceed, a continuum description is introduced, in the form of equations of motion for the coarse-grained hole and 
charge density profiles. Since the latter are both conserved quantities, the equations of motion take the form of 
continuity equations. They can be derived phenomenologically [ p5| or directly from the master equation [ p^ . In the 
latter case, we first write a set of equations for the local averages, {4>x,y) and {ipx,y) on discrete space (with lattice 
constant 1) and then take a naive continuum limit, e.g., we approximate (^^{(f>x+i,y — 4>x-i,y)) by a first derivative 
with respect to x, etc. A mean-field assumption is necessary since two-point correlations must be truncated in order 
to obtain a closed set of equations. These can easily be written in general dimension d: 

9t(/.(f,t) - V - {v0 + £0Vy} 

dtW,t)^^ ■{<t>^^-i^^4>-£4>{i~4>)v], (9) 

Here, the hole density (j) and the charge density ijj are functions of the d-dimensional coordinate r (with associated 
gradient V), and time t. The drive appears in these equations via its coarse-grained equivalent, the effective drive £: 

£{E) = 2 tanh(£;/2). (10) 

pointing along unit vector y. A diffusion constant has been absorbed into the time scale. Derivatives higher than 
second order have been neglected, anticipating smoothly varying solutions. The equations have to be supplemented 
with periodic boundary conditions and the constraints on total mass and charge. For later reference, we also define 
the parameter 

e ^ £Ly (11) 

which will play the role of a scaling variable. 

Time-independent solutions of these equations reflect stationary phases of the discrete model. The disordered 
phase corresponds to a homogeneous solution, which is stable with respect to small perturbations provided m does 
not exceed a threshold value m/f , given by 

mH = [l + {2TT/ef]/2 . (12) 

The profiles in this phase are uniform. In our case, where the lattice is nearly completely filled, i.e., m < 1, we need 
e < 27r in order to find a stable homogeneous steady state. For an electric field E = 1.0, this implies rather small 
system sizes {Ly < 7). 

To find a steady-state solution which corresponds to a transverse strip, we seek solutions that are inhomogeneous 
in the y-coordinate only. Eqns (||) can be integrated once, with integration constants being the hole and the charge 
currents. The former vanishes by symmetry at zero total charge. The latter, being non-zero in general, will be 
denoted by j£. After expressing -ip in terms of (j): 

and rescaling the spatial variable to z = y/Ly, we obtain an ordinary differential equation, for the function x= ^/4'- 

x"{z)/e^^~jxHz)+x{z)-l. (14) 

To satisfy the boundary conditions, x should be periodic with period 1. Writing x in terms of a potential (l/e^)x" = 
— {d/dx)V{x), a further integration leads to x' = e\/2(C/ — V), where U is another integration constant. Unique 
solutions exist for j < 1/4 and appropriate U. Introducing the three roots xi !i X- ^ X+i defined via 2(U — V{x)) = 
(2j/3) {x+ — x){x~ X-){x~ Xi)i the solution [f^ can be written using Jacobian elliptic functions [p7|: 



4 



X{z) = X+ - (X+ - X-) sn^ (e^V 076) (x+ - Xi)) ■ (15) 

in the interval < z < 1/2. The other half of the interval, 1/2 < z < 1, is described by symmetry around the point 
z = 1/2. Thus, the hole density takes its minimum at (f>{0) = 1/x+j ^ind its maximum at 4>{^) = 1/x-- The third 
root, xii li6S outside the physical domain. It is convenient to define the parameters p and R: 

P= (x+-X-)/(x+-Xi) (16) 
R^i'iKip)/e)\ (17) 

Here, K stands for the complete elliptic integral and is a function of p. Quantities of interest, such as the mass 
m or current j, can be expressed in terms of p and R: 

l-4j=i?2(l-p + p2) (18) 

, _ [l-R^{l-p + p^)] n{n\p) 

2(1 +R + pR) Kip) ' ^'^^ 

where n(n|p) is the complete elliptic integral of the third kind and n = 3pi?/(l + + pi?). In principle, Eqns ( p^ ) and 
( |l9| ) can be inverted to give the physical parameters m and e in terms of R and p. In practice, it is easier to generate 
functions of interest, e.g. j(e, m) or the order parameter Ql{£, "m), parametrically in p, which is allowed to range from 
to an upper limit po{e) < 1 The upper limit po is defined by the vanishing of the current, j{e,Po) = 0, and 

plays a particularly important role in the context of this study: According to Eqn (p^, the mass density m tends to 
unity as p approaches its upper limit po(e)- Thus, only values of p near pq will be of interest here, since our focus is 
on nearly filled systems. This observation is used later for approximations. 

The solution for xiz) generates both, hole and charge, densities: cj) = l/x(^) a-nd ip = x'/(^x)- These solutions 
describe the ordered phase, i.e., particle-rich strips transverse to the field. For fixed mass, they depend only on the 
parameter e = £ Ly and the variable z — y/Ly; thus, these functions satisfy scaling in these variables. Moreover, the 
order parameter Q^, is a function of e alone, since the spatial variable is integrated out. Here, however, we have to 
be rather more careful: since our system, irrespective of its size, will always contain only a single hole, the mass is 
inherently size-dependent. We will return to this issue in the next section. 

Since it is cumbersome to work with Eqn ( p^ ) directly, its approximation for e ^ 1 is very useful p^ . The 
sn-function can be replaced by a tanh-function, and the argument simplifies : 

X(z)~x+-(X+-X-) tanh2(ez/2). (20) 

As a result, a (weak) discontinuity appears in the first derivative of x at the symmetry point z — 1/2. This is unfor- 
tunate for our purposes, since z — 1/2 is also the location of the maximum hole density. A different approximation, 
to be presented in the next section, resolves this difficulty. We note in passing that Eqn (^o|) takes the form of the 
soliton in the Korteweg-de Vries equation pst . 

Clearly, one should not expect such a mean-field theory to provide a quantitatively correct description of the phase 
transition. However, it gives excellent qualitative insight into the instability and the phase diagram JlTf . Moreover, 
since our interest here focuses on behavior deeply in the ordered phase, fiuctuations do not play a significant role, and 
a mean-field theory should be very reliable. In fact, we will see that its main limitations do not arise from the neglect 
of correlations, but from taking a naive continuum limit. 

We conclude this section with a few technical details of the simulations. The linear system sizes, and Ly, range 
from 16 to 48, with E in the range 0.2 to 1.2. A characteristic parameter set is that of our "reference system", which 
will appear in all scaling plots: E = 0.8 and x Ly = 16 x 24. Thus, the mass density differs from unity by at most 
0.4%, so that m = 1 is often an excellent approximation. The statistical error of the simulation results is of the order 
of 5% and thus much larger. All initial configurations are random. In one MCS, a nearest neighbor of the vacancy 
is chosen at random and an exchange is attempted with the rates (Q). Averages are computed from 100 independent 
samples for each choice of parameters. The approach to steady state is extremely slow |2^ for larger system sizes and 
sets real-time limits on our simulations. For example, a system with 16 x 24 sites at E = 0.8 requires approximately 
5 X 10^ MCS to reach the steady state. If we increase Ly from 24 to 36, which is a factor of 1.5 in system size, the 
required number of MCS increases by roughly a factor of 10. 

While averaging, e.g., Ql, is rather simple, by first measuring Ql for each sample and then averaging these data, 
some effort is needed to compute average density profiles from the configurational data. Due to translational invariance, 
strips can be centered at any y, and a careless average would "wash out" any inhomogeneities. To avoid this, we first 
shift the ordered strips in the different samples in such a way that they match before we average. A natural choice 
would be to center all strips on, e.g., y = 0, by normalizing the phase of the largest wavelength Fourier component 
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of the profile | |19| . This is particularly useful when profiles are measured near the phase transition. Here, however, 
we will mostly take data deeply in the ordered phase, where the hole is essentially trapped. Thus, for each sample, 
we keep track of the j/-position of the hole and determine the maximum of the hole density after a large number 
of MCS. This maximum marks the interface between the positively and negatively dominated regions, the former 
located "up- field" from the latter. The charge density profiles from different samples are now shifted such that these 
maxima coincide, and averages can be taken. Clearly, this procedure would run into difficulties if the interface were to 
wander significantly while the data for the hole density profile are being accumulated. However, for the choices of the 
control parameters considered here, this does not appear to present major problems since fluctuations of the interface 
position are rather small. Moreover, they are very slow; thus, the time scales over which the interface remains well 
localized are sufficiently large to determine the maximum of the hole density very precisely. 



III. SCALING BEHAVIOR IN THE STEADY STATE 



As an introduction to the discussion of scaling properties, we illustrate the process by which the system approaches 
the steady state. A series of snap shots, taken at different MC times, demonstrates, first, why the dynamics is so slow, 
and second, already suggests one of the key hypotheses of this work, namely, that the steady-state interfaces are well 
separated from one another. Figure 1 shows this series for our reference system. The negative (positive) particles are 
colored black (white) and the empty site is marked gray. A coordinate system is introduced in the usual way, i.e., the 
x-direction lies horizontal, the y-direction vertical and the i?-field points upwards. 

Starting from a random configuration (Fig. la), the system remains disordered for early times (Fig. lb). Eventually, 
by allowing positive (negative) particles to move preferentially upwards (downwards), the hole begins to segregate the 
two species. The early stages of this process are discernible in Fig. Ic, where an interface between regions of opposite 
charge begins to develop. The position of this interface is determined by random fluctuations in the system. Clearly, 




d) e) f) 

FIG. 1. Snapshots of an Lx x Ly system witli = 16 and Ly — 24: a.t E — 0.8, at different numbers of MCS: a)l, b)10'^, c)10'', 
d)10*, e)10^, f)10^. The initially disordered system undergoes a charge segregation. In the ordered steady-state configuration 
the two oppositely charged particles are separated by two different interfaces. Minus particles are colored black, plus white and 
the hole is gray. i5-field and y-direction point upwards. 
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due to the periodic boundary conditions, a second interface must also form. After 10^ MCS (Fig. Id) the segregation 
of charges, and hence the two interfaces, are quite apparent. Due to the drive, the hole moves rapidly to the top 
(bottom) in regions of predominantly negative (positive) charge. Thus, it tends to remain near the interface which 
separates positive particles on the top from negative ones on the bottom (the lower interface in Figs. Id-f). In 
contrast, it is rapidly driven away from the opposite interface. We will refer to the former (latter) interface as the 
"downstream" ("upstream") one. The two interfaces are well separated, for this choice of parameters, and exhibit 
rather different morphologies: The downstream interface is quite sharp, while the upstream interface appears to be 
much more diffuse. To increase the degree of order in the system, the hole has to travel to the upstream interface 
before it can move another charge to a preferred position. Since this requires a series of field-suppressed jumps, the 
approach to the final steady state is very slow. As the ordered domain surrounding the downstream interface grows, 
the hole becomes strongly localized. A quantitative analysis of this ordering process will be provided elsewhere |22j . 

A picture of a typical steady-state configuration is shown in Fig. If. To characterize these structures, we investigate 
three characteristic quantities: the order parameter , which provides a global measure of order, as well as the average 
hole and charge density profiles which carry more detailed information about ordered configurations. All three of them 
are easily computed within the mean-field theory, as we shall presently see. Since the spatial inhomogeneities form 
along the y-direction, the system size is not expected to play an important role. Simulations confirm this, provided 
the aspect ratio L^/Ly does not exceed a certain threshold value which is at least 6 in our case. For larger aspect 
ratios, strip configurations with nonzero winding number may begin to form, introducing an Laj-dependence into the 
problem [|6|. These, however, are not the subject of the present study. 



A. The Order Parameter 



We begin by calculating the order parameter in the mean-field approximation. Starting with the definition (|g), 
we first express it within the continuum theory. Clearly, the summation over the sites in the y-direction should be 
replaced by an integration. Using the rescaled variable z = y/Ly and exploiting the symmetry of the profiles around 



z = i, we obtain: 



!L = — 



1/2 2 /-i/^ v'2 

mz)f dz^— ^dz . (21) 

In the last equality, we have recast ip in terms of x- To proceed, we change integration variable, from z to x- The 
limits of the integral are transformed to x(0) = x+ and xi'^) = where we recall that these are zeroes oiU — V{x)- 
With x' = V(2je2/3) (x+ - x){x - X-)(x - Xi), we find 

4j 



, (Qi + Q2 + Q3 + Qi) (22) 
dm 



where we have introduced 



Qi = - ^dx (23) 
Jx- X 

Q2 = (xi X-X+) -^dx (24) 
Jx- XX 

rx+ I 

Qs = iXi + X- + X+) / —dx (25) 
Jx- X 

Q4^-{X+X-+X+Xi+X-Xi) — jdx . (26) 

Jx- XX 

The last two integrals are evaluated easily, giving 1/2 and (1 — m)/2. The first two integrals can be reduced to 
complete elliptical integrals of the first, second and third kind. 

The resulting expressions can be expressed in more compact form, using the parameters rn, p and R (see Eqns ( [l^ ) 
and ([l^)). For that purpose, it is helpful to write the three roots in terms of p and R: 

_ 2{l + R{p-2)) 

l~R^l-p + p^) ^ 
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_ 2(l + i?(l-2p)) 

1 — (1 — p 



Then, we invoke Eqn (18) to replace the current j and Eqn ( pl^ to eliminate the elliptic integral of the third kind. 
Collecting, we obtain Q^: 

1 .. . 



0. = l--|i?(p-2) + l + 3i^^|. (30) 

E and K are the complete elliptical integrals of the first and second kind According to Eqns ( p^ and ([l9|), Ql 
is a function of p and e only which can be generated parametrically in p. 

So far, our discussion is valid for arbitrary mass m. Let us now consider the case of a single vacancy, namely 
m = 1 — l/{LxLy) which corresponds to a system near complete filling. Considering only the leading terms in 
an expansion in powers oi 5 = l/{LxLy), the left hand side of Eqn ( |l9[ ) is just S. As a consequence, the factor 
1 — R^{1 — p + p'^) on the right hand side is 0{5), and so is the current j, given by Eqn (p^. Recalling that the 
upper hmit po of the p range is defined by the condition: j(e,pQ(e)) = 0, we conclude that p = po + 0{6) for our case. 
Tracking the effect of the finite-size corrections through our preceding calculations, we find that xii X+ ^-iid X- s-rs 
all 0{LxLy), by virtue of their common denominator. To leading order, the hole density (j) = 1/x is therefore 0{S) 
as one might have anticipated! In contrast, the charge density is 0{1), due to Eqn (|l^). Since 6 is very small in our 
study, all but the leading terms will be neglected in the following. Then, becomes a function of e alonel This 
prediction is easily checked by Monte Carlo simulations. 

In this spirit, we invoke Eqn ( p^ for j = and rewrite it as — [AK {po)]'^ {1 — po +_Po)- Now, e can be computed 
numerically for a set of discrete values oi p in the interval [0, 1]. The values of (5i(e, m = 1), derived in this way, are 
shown as the solid theoretical curve in Fig. 2. 

For large e 3> 1, the approximations R ~ \/l — Aj and p ~ 1 are valid. Within the same approximation, j can be 
replaced by j = exp(— 77ie/2) which is vanishingly small. With to = 1, simplifies to 

Ql = 1 - ^ (31) 

This gives rise to the dashed curve in Fig. 2. Comparing this approximation to the exact mean-field result, we see 
that both expressions are indistinguishable for e > 18. For smaller values of e, the approximation underestimates the 
order parameter slightly. 

Turning to simulation results, we first test the expected scaling in e. Figure 2 shows data for the order parameter 
Ql, for different square systems (ranging from 20 x 20 to 35 x 35) and different electric fields [E = 0.2 to 1.0), plotted 
versus e. Each data point is an average over 30 - 50 samples. The size of the error bars is about 0.05 units. 



T3 




20x20 
25x25 
30x30 
35x35 
1-6/e 
QL(mft) 



30 



8=eL„ 



FIG. 2. Plot of order parameter vs e for different square system square sizes (20 < Lx = Ly < 35) and different electric fields 
(0.2 < E < 1.0). The upper line shows the result of from mean field theory, the lower line is its approximation 1 — 6/e. 



8 



Within the accuracy of our data, all points lie on the same curve, corresponding to the scaling function. The latter 
appears to tend towards zero for e < 6. This is consistent with the stability limit of the inhomogeneous solutions, 
Eqn (p^), which implies that for m ~ 1, an inhomogeneous solution can exist only if e > 2tt. Once the transition to 
the homogeneous phase has occurred, the order parameter is of the order oi l/L^ < 0.05. With increasing e the order 
parameter approaches its upper limit, i.e., 1. 

It is remarkable, however, that the large e approximation produces a better fit to the simulation data than the exact 
mean-field result, especially for e < 10. It is conceivable that the intrinsic errors of the mean-field theory approach 
are partially compensated by the large e limit. Further studies are required to test this possibility. Focusing on the 
region e > 15, the simulation results all lie about 0.02 units above the theoretical curve. While these deviations are 
within the error bars of the data, they are too systematic to be ignored. Closer scrutiny reveals that the results of the 
large system sizes tend to be closer to the theoretical curves than those for small system sizes, which indicates that 
the differences between simulation and mean-field results are at least partly due to finite size effects. We will return 
to this question at the end of Section III.C. 



B. Charge and Hole Density Profiles 



While the order parameter carries only global information about spatial inhomogeneities in the system, the charge 
and hole density profiles retain far more detail, allowing us to distinguish the oppositely charged domains and their 
interfaces. Based on the mean-field theory, we expect these densities to satisfy scaling in e and z. This is borne out 
by the simulation results which are presented in this section. 

In order to exhibit the scaling of the densities, four different parameter sets {E, Ly) are simulated, generating 100 
samples for each. The system length Ly and the electric field E are varied in such a way as to keep the parameter 
e constant at 18.24. To avoid unwanted cancellations, we shift the maximum of the hole density in each run to 
z — Q before averaging. The charge profiles are shifted accordingly. Thus, z covers the interval (—0.5,0.5), and the 
"downstream" interface is centered at the origin. In addition, we normalize the hole profile in such a way that all 
profiles enclose the same area. 

A comment on this normalization is in order. Recalling the constraint on the total density, we have I ~ m — 
l/{LxLy) = (l/Ly) J^^ 4>{y)dy for a single hole. Thus, we have I = J^^ 4'iy)dy = L^Ly (j){z)dz so that Lx(f>{y) 
can be interpreted as the probability density for finding the hole in row y. Similarly, LxLy(f>{z) is the probability 
density for finding the hole at position z. Thus, normalized plots for the hole density show the associated probability 
density, and the area under each curve is just 1. Moreover, since (/^{z) = 0{l/{LxLy)), according to the finite-size 
analysis in Section IIIA, the normalized quantity depends on z and e alone. No such normalization is required for the 
charge density, since it is already of 0(1) in the system size. 

To test for the anticipated scaling in e and z, Monte Carlo data for the (normalized) hole density profile are 
presented in Fig. 3a, and the charge density profile is shown in Fig. 3b. Since all data points collapse onto the same 
characteristic scaling curve for hole and charge profiles, respectively, the theoretical prediction is clearly confirmed. 
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FIG. 3. Scaling plot of the hole (a) and charge (b) densities, 
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Beyond demonstrating scaling, these plots provide a more quantitative characterization of the spatial structures in 
the system. Since e here is the same as in Figs, la-f, Fig. 3b shows the associated steady-state charge density profiles. 
We can see clearly that the particles are ordered in two regions, filling the whole system. Each of these regions consists 
essentially of one species. They are separated by two interfaces. The maximum of the hole density lies at the center 
{z = 0) of the much sharper downstream interface where the field tends to localize the hole, while the minimum of 
the hole density marks the more diffuse upstream interface. 

To explore the size- and field-dependence of our system further, it is interesting to vary the system length Ly and 
the electric field E independently, not keeping e constant. Of course, given the excellent data collapse of Figs. 3 and 
4, we cannot expect global scaling over the whole y-range. We will see, however, that certain regions of the profiles, 
centered on the two interfaces, still scale. 

We first report simulations at constant electric field E = 0.8 and transverse size — 16, increasing the longitudinal 
system size Ly from 20 to 32 in steps of 4. The (normalized) hole densities observed in these simulations are summarized 
in Fig. 4a, plotted vs y rather than z ~ y/Ly. 




We observe that the graphs associated with different Ly span different ranges of y, but are otherwise essentially 
indistinguishable in the central region. Thus, the width and the maximum of the hole density, and hence the width 
of the "downstream" interface, are not affected by changes in the longitudinal system size when plotted vs the y- 
variable. We conclude therefore that the characteristics of this interface are controlled by the electric field alone (cf. 
next section). 

This behavior is also borne out by the charge density profiles. Fig. 4b. According to Eqn (p^), the steady-state charge 
and hole density profiles are related via ip{y) — —(p'{y)/{£(j){y)). Thus, the charge densities near the "downstream" 
interface should also be independent of Ly, in agreement with Fig. 4b. On the other hand, the regions of nearly 
constant charge density must broaden, to reflect the increasing system size. Thus, the profiles do not collapse at 
the edges of the plot. However, the similarity of their form near y = ±Ly/2 suggests that the upstream interface 
might scale also, provided the profiles are shifted appropriately. This is indeed confirmed by the simulations (cf. 
Fig. 5b below). Thus, the slopes and widths of the profiles, near both interfaces, are independent of system size. 
The remaining effect of Ly is very simple and can be observed in Fig. 4b: Outside the interfacial regions, the charge 
densities saturate very rapidly at ±1, and these saturated regions expand or contract to accommodate the selected 
system size. 

It is now quite apparent how the profiles should scale if Ly remains fixed and E is varied instead. Since the 
interfacial regions are independent of Ly, but scale in z and e, they must depend on y through the combination £ y. 
To check this conjecture, we fix the system size at L^x Ly = 16x 24, while the electric field increases from 0.4 to 1.2 in 
steps of 0.2. Fig. 5a, shows the (normalized) hole density profile plotted vs the scaling variable £y. The data collapse 
in the interfacial region is excellent, except for the smallest field E — 0.4. This value of E, however, is rather close 
to the transition line where the mean-field theory is likely to break down. Focusing on the larger £"s, it is apparent 
that the width of the downstream interface scales as l/£. Turning to the charge densities, Eqn ([l^) implies that ip{y) 
is also a function of £y near the downstream interface. This is indeed confirmed by simulations. To illustrate the 
scaling of the upstream interface, we present Fig. 5b: Here, all profiles have been shifted by Ly/2, in order to center 
the upstream interface at the origin. Clearly, this interface also scales in the variable £ y. For completeness, we note 
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that the only profile that does not reach saturation, is the one for the smallest E, since this value is quite close to the 
phase transition. 

Let us summarize the key findings of the simulations. First, the data for each profile collapse onto a single, global 
scaling curve if plotted as a function of z = y/Ly at constant e. Moreover, focusing only on the interfacial (as opposed 
to the saturation) region, we find that both interfaces are independent of Ly and that their widths scale as 
provided we are not too close to the transition to the homogeneous phase. In the next section, we will consider these 
findings in the light of our mean-field theory. 
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FIG. 5. Plot of the hole (a) and charge (b) densities for a range of E, vs £y. 
to the minimum of the hole density (upstream interface). 



16. Note that in (b) y = corresponds 



C. Independent Interface Approximation 

In the following, we present an analytical description of the interfaces which is then tested by detailed Monte Carlo 
simulations. We will focus in particular on the charge density since it directly determines the order parameter. First, 
we invoke the "large e" approximation to describe the upstream interface. A complementary approximation, 
involving a different version of the large e limit, is then developed to describe the other (downstream) interface. The 
key assumption here is that the two interfaces are independent from one another, i.e., they are separated by sufficiently 
large "saturation" regions: regions which are entirely filled by either positive or negative charges, so that ip = ±1 
there. Such profiles result provided e > 18, as demonstrated by Figs. 3b, 4b, and 5b. The associated hole densities 
are approximately zero except near the downstream [z — 0) interface (cf. Figs. 3a, 4a and 5a). 

Returning to the large e limit of our mean-field theory, we recall that e ^ 1 is equivalent to p ^ 1. In this limit, we 
may replace the Jacobian elliptic function sn by tanh |17|. In practice, this is already a good approximation for e > 15. 
It is easy to check that this results in (mean-field) profiles with vanishing hole densities near the upstream interface 
and saturated charge densities between the interfacial regions. Thus, this limit is consistent with our assumption of 
"independent" interfaces. 

To describe the upstream interface, we start from Eqn ( |20| ) for large e: 

X(z) = x+ - (X+ - X- ) tanh' (ez/2) (32) 

Note that, due to the symmetry of x, this equation holds for the interval (—1/2, 1/2). At z = 0, x takes its maximum, 
so that this is a good approximation for the upstream interface where = 1/x is minimal. The largest deviation from 
the exact mean-field solution occurs at the boundaries, i.e., near the downstream interface, since this approximation 
violates the periodic boundary conditions: x'(~l/2) ^ x'(l/2)- The current is exponentially suppressed for large e, i.e., 
j ~ 66"™*^/', and x+ and X- can be expressed in terms of j, namely, x+ — X- — JjV^ ~ 4j and x+ — ^(l-l-2-y/l — 4j) 
|p7| . The hole and charge densities are now easily derived. In particular, recalling that y = LyZ, we can already 
read off the width ^„ of the upstream interfacial region: ^„ — 2Ly/e ~ 2/£ which is consistent with the data. More 
specifically, we can compute the charge density from Eqns ( p^ and (|3^). Neglecting terms of 0{j), we find 

^-(2/) = -tanh(£j;/2). (33) 
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At the boundaries, Eqn ( |3^ ) results in ■0(±Lj,/2) =f1 in the large e limit, which confirms that this approximation 
violates the boundary conditions. However, it does describe the interfacial region near y = very well. Rather than 
quoting the hole density explicitly, we only note that it is very small near the origin, namely 0{j). 

In order to capture the downstream interface, we introduce another method. Since j = to excellent accuracy 
for large e [ p^t , we return to the mean-field equations, (j^), and integrate them, setting both integration constants, 
i.e., hole and charge current, to zero. This is actually an equilibrium approximation, as we shall discuss below. The 
simulation data suggest the boundary conditions: (j){±Ly/2) ~ and ip{^Ly/2) ~ ±1. In this approximation, the 
downstream interface, corresponding to the maximum of the hole density, is localized at the origin. Written in terms 
of the variable y, Eqn. ( p^ ) for x simplifies to 

x"iy)/£' = x{y) - i- (34) 

This is easily solved, subject to the specified boundary conditions: 

X(y) = 1 + ccosh(£y) . (35) 
To ensure that the hole density is strictly positive, we demand c > 0. This constant can be determined explicitly from 

Ly/2 

'iy/2' 



the mass constraint, namely \ = Lx /^2^/2 4'{v) dy, whence 



^LxL^ arccos ( - ) = 1 (36) 

in the large e limit. In our simulations, is at least 16 and £ at most 2, so that c > 24 follows. Thus, we can expand 
Eqn for large c, resulting in 



LrL 

-IT 



(37) 



In fact, this approximation is already very good for c > A. Next, we compute the charge density, 

'yf^^„i ^tanh(gy) . (38) 

The last approximation is very accurate since c > 24. Again, we can read off the width of the interfacial region, 
= \/£. Similar to the downstream interface, the width scales with 1/f , in agreement with the data. Intriguingly, 
however, our approximation is capable of reproducing the observation that the downstream interface is narrower than 
the upstream one. Whether the measured widths differ by a simple factor of 2, as predicted by our calculation, awaits 
a more quantitative comparison with Monte Carlo data. 
In contrast to the upstream interface, (j) is nontrivial here: 

'^^^^ " l + ccosh(£2/) ~ c cosh(£y) ^^^^ 

confirming the width of the downstream interface = 1/f • Away from the origin, the hole density again decays very 
rapidly, to match with its value near the upstream interface. 

Before turning our focus on a comparison of these results with computer simulations, a last remark on the ap- 
proximation of the downstream interface is in order. Imposing brick- wall (i.e., closed) rather than periodic boundary 
conditions, the approximation taken here (setting the current to zero) becomes exact. Moreover, the brick-wall sys- 
tem is an equilibrium one. The hole will accumulate positive (negative) charged particles at the top (bottom) of the 
system, thus establishing our boundary condition ip{±Ly/2) —^ ±1. Clearly, only one nontrivial interface remains in 
this case, namely, the downstream one. In the steady state, the bias traps the hole near this interface. This fixes the 
boundary condition (j){±Ly/2) —^ 0. In this sense, our approximation for the downstream interface is equilibrium-like. 



Returning to our model, we have obtained two compact equations, (33) and (38), for the charge density. Since 
ijj^y) ~ ±1 between the interfaces, to excellent accuracy, the whole system can be described in terms of the two 
interfaces, provided we match them appropriately. As an example. Fig. 6b shows a 16 x 24 system with E = 0.8. 
The data points result from a Monte Carlo simulation while the solid and dashed lines refiect our two interface 
approximations, Eqns ( ^3| ) and (^), respectively. For the narrower interface (downstream; in the center), the match 
is nearly perfect, while for the wider interface (upstream; at the edges of the figure) the slope of the tanh-function is 
slightly too small compared to the computer simulation. The agreement is nevertheless remarkable. 
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24 and E — 0.8. The dashed hne denotes the approximation 



FIG. 6. (a) Plot of the hole density (+) for = 16, 1 
Eqn (^^, for the downstream interface, for c = 65.58. Note that the hole density is vanishingly small near the upstream 
interface. (b) Plot of the charge density (+) for Lx ~ 16, Ly — 24 and E — 0.8. The dashed and solid lines are the two interface 
approximations, Eqns (p3) and (pq), matched at y — ±4. 



Next, we compare the approximation for the hole density with MC resuhs. Here, we use Eqn p9| ) for the whole 
system since the hole density is vanishingly small except in the central region of the downstream interface. Simulation 
data and the analytic approximation for cf) are presented in Fig. 6a. A small quantitative discrepancy is observed at 
y = OjLe., the center of the downstream interface, whereas all other data points are remarkably well reproduced by 
Eqn®. 

Given the results for the interfaces, we finally return to the order parameter Ql. Here, we will see that the 
independent interface approximation provides us with a very intuitive picture for the approximate form (31). Since 
the steady state exhibits complete order in one region of positive and another one of negative particles, the deviation 
of Ql from unity originates near the interfaces. We can easily compute the contribution to Ql for each interface 
separately, using Eqn (pTJ). The wider (upstream) interface reduces Ql by A/{£Ly), while the narrower downstream 



interface lowers it by ), resulting in a net Ql = 1 — 6/{£Ly), in agreement with Eqn (^l]). Thus, this for 

simply tallies up the contributions of two well-separated interfaces, while the fully saturated regions give rise to the 
1. 

While these data for the scaling of order parameter and profiles are very convincing, the question of their range of 
validity must be raised. First, we should anticipate a breakdown of mean-field theory near the onset of the transition 
to the uniform state. This limits our analysis to e > 6, corresponding to, e.g., E > 0.25 for a system with Ly = 24. 
For larger values of e (but below an upper limit to be discussed shortly), scaling in e and z is observed to hold. In 
order to have well-separated interfaces, we also require e > 15. Beyond this threshold, the interfacial regions of the 
profiles scale very cleanly with l/£. 

In addition to a lower limit, there is also an upper limit for our analysis. Recalling Eqn ([To|), the effective drive 
is bounded: £ < 2, due to the tanh-function, even for very large values of the microscopic E. Thus, within our 
mean-field theory, the interfacial widths cannot become arbitrarily small. For example, for the narrower downstream 
interface ■Cd < 2 in units of the lattice spacing, and i^{±2) — ±0.96, from Eqn (^). Thus, mean-field profiles require 
at least 4 lattice spacings, to interpolate between the fully saturated regions. In contrast, measured charge density 
profiles for large E (e.g., E — 2.0), jump from —1 to -1-1 over just two lattice spacings! Such profiles are so sharp 
that our continuum limit fails to reproduce them: they can hardly be considered smoothly varying functions. As a 
result, the mean-field theory underestimates the order parameter for large values of E which explains the systematic 
deviations of the smaller system sizes in Fig. 2. For example, e = 20 in a 20 x 20 system corresponds to E — 1.2 
where this phenomenon is already noticeable. At a purely phenomenological level, we can extend the validity of our 
mean- field description if we retain the form of our equations, (|^) , but replace the effective drive £ by the microscopic 
field E everywhere. Mathematically, this requires keeping explicit track of the lattice constant a, followed by taking 
the hydrodynamic limit [p9[|, i.e., a — *■ at fixed drive, system size and mass. Since the lattice constant appears in the 
rates, Eqn (|l]), the effective drive takes the form £ — 2 tanh(i?a/2). In the original discrete version of Eqn (H), the 
lattice constant appears in terms such as £ (^\{(j)x+a.y — 'Px-a.y)}- In the limit of vanishing a, this expression simplifies 
to Ea^d(t>/ dx. Since the diffusive terms (e.g., V^0) also generate a factor a^, the latter can be absorbed into the time 
scale so that we recover Eqn (O), with £ replaced by E. Thus, all of our analytic results carry over, provided E takes 
the place of £ everywhere With this modification, the agreement of MC data and analytic description extends 
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to the largest fields studied, namely, E = 2.0. 

To some extent, even the measured profiles do not reproduce the actual sharpness of the data fully. Since the 
downstream interface can form at arbitrary location within the lattice, one should allow for noninteger shifts, i.e., 
shifts between and 1 (a = 1) modulo multiples of the lattice spacing, in order to produce accurate averaged data. 
This subtlety is not accounted for in our simulations, as seen from the discussion at the end of Section Thus, 
the actual interface is slightly smeared out when we average profiles by superposing the maxima of the hole density. 
Details can be found in Ref. To summarize briefly, our mean- field theory, in the form of Eqns (||), gives excellent 
results provided e ^ 6 and E < 1.0. If a systematic hydrodynamic limit is considered, the validity extends further, at 
least to E < 2.0. 



IV. CONCLUSIONS 

In this work, we focused on the scaling behavior of ordered steady states in a simple lattice model. A fully periodic 
lattice is filled with equal numbers of positive and negative "charges" , except a single site which remains empty. An 
external "electric" field, applied along one of the lattice axes, biases the motion of the particles. The dynamics is 
vacancy-mediated in that only vacancy-charge exchanges are allowed. The particles interact only through an excluded 
volume constraint. 

This system develops spatial structures if ELy, i.e., the product of drive and system size along the field direction, 
exceeds a critical value. Then, a charge segregated strip, oriented transverse to the field, forms around the hole and 
grows until it fills the whole system. The two oppositely charged regions are separated by two interfaces with distinct 
characteristics: One interface, the "downstream" one, attracts the hole, the other (upstream) repels it strongly. This 
asymmetry finds its origin in the charge separation induced by the external field: while the hole moves rapidly along 
the field in the negative region, its preferred direction is reversed in the positive region. 

Continuing earlier studies |l5|Jl7[] , we investigate the scaling properties of an appropriate order parameter and the 
hole and charge densities, as the external control parameters E and Ly vary. The transverse system size plays 
no role except in finite-size corrections. Monte Carlo data are compared to the predictions of a mean-field theory 
in which the drive appears through the effective parameter £ = 2tanh(_E/2). The agreement is excellent, provided 
e = £Ly ^ 6 so that we are in the ordered phase, and E < 1.0 to maintain fairly smooth profiles. The transverse 
system size Lx plays no role except in finite-size corrections. In particular, we can describe the charge density profiles, 
with remarkable accuracy, in terms of two non-interacting interfaces, separated by perfectly ordered regions. The 
interfaces themselves are determined by the drive alone, independent of system size, and their widths scale with l/£. 
For fields i? > 1.0, the data show very steep slopes in the interfacial regions which cannot be captured correctly by a 
naive continuum limit. Remarkably, the mismatch between data and mean-field theory is significantly reduced if we 
substitute the microscopic field E for the effective £ in the (mean-field) interface approximations. The emergence of 
the latter can be understood in the limit of vanishing lattice constant. We should caution, however, that this limit 
must also eventuelly break down since it does not commute with the limit E oo. Since the details of the continuum 
limit appear to play a key role here, it would be interesting to analyze the discrete precursor of Eqn (^). In this case, 
the natural parameter should be £ alone. 

Another interesting question concerns the character of the interfaces when a finite density of vacancies is present. 
In this case, the downstream interface "splits" into two halves, separated by an empty region. Clearly, in addition to 
E and Ly, the overall mass m now enters the criterion for having independent interfaces. Provided the appropriate 
condition is met, however, we expect that the interfacial profiles still depend only on £: local structures appear to be 
controlled entirely by the drive. 

Finally, our study prepares the ground for the exploration of dynamic phenomena in driven two-species models. 
Having established the scaling properties of the final steady states, work is in progress to investigate how they develop 
from random initial conditions . 
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